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L  Introduction 


From  February  -  March,  1992  Professor  Ernie  Tuck  from  the  Applied  Mathematics  Department  of 
Adelaide  University  undertook  an  investigation  for  DSTO  of  Submarine  Wakes,  The  emphasis  of 
the  investigation  was  to  compute  the  internal  waves  generated  by  the  passage  of  the  submarine 
through  a  stratified  ocean.  Ebor  Computing  was  engaged  to  provide  software  support  to  Professor 
Ernie  Tuck  and  DSTO.  This  support  consisted  of  developing  software  to  model  these  internal 
waves  and  produce  output  in  a  form  so  that  it  may  be  srudied  or  used  by  graphics  packages  such  as 
MATLAB  and  PLOTZ  to  display  images  of  the  wakes. 


Figure  1  is  a  diagram  showing  the  overall  structure  of  the  software  developed. 


Figure  1  :  System  Diagram 


This  document  describes  the  program  developed  by  Professor  Ernie  Tuck.  Dr.  Graham  Fumell,  Mr. 
Tony  Legg  from  DSTO  and  Mr.  Michael  Carroll  of  Ebor  Computing. 


The  program  was  developed  on  a  SparcStation  using  MATLAB  and  PLOTZ  to  generate  plots  of  the 
wakes.  To  run  the  program  you  need  to  copy  all  the  source,  parameter  and  density  p  ofile  files 
into  your  directory  and  compile  them,  see  Chapter  5  on  Compiling.  You  will  also  need  access  to 
the  IMSL  Library  Routines  (Fortran  Version).  Once  the  program  has  been  compiled  the  program  is 
run  as  follows  : 

>  wake 

2.1.  Parameter  File 

The  user  is  not  prompted  for  any  inputs,  all  inputs  to  the  program  are  contained  in  the  parameter 
file  VVAKE.PAR.  Figure  2  shows  a  sample  parameter  file  : 

NB  :  The  Parameter  file  must  have  the  following  format  : 

3  Header/Comment  lines  :  Gives  information  such  as  file  purpose  and  date. 

23  Parameter  Lines  :  First  71  Characters  contain  a  comment  on  the 

parameter,  the  rest  of  line  contains  the  actual 
parameter. 


This  file  is  read  by  WAKE  and  contains  initial  values  for  some  of  WAK 

S '  s 

options.  Michael  Carroll  13.4 

92 

NUMTHETA  -  The  number  of  theta  values  to  calculate  data  for  (<=30 ) 

30 

MODE  -  The  mode  to  calculate  data  for 

1 

THETAMAX  -  The  maximum  theta  value  allowed  (<=  Pi/2) 

1.56 

SUBDTH  -  The  depth  of  the  submarine  (in  metres) 

50. 

SUBSPD  -  The  speed  of  the  submarine  (in  metres /second) 

2.5 

SUBLEN  -  The  length  of  the  submarine  (in  metres) 

100. 

SUBRAD  -  The  radius  of  the  submarine  (in  metres) 

5  . 

XSTART  -  The  starting  point  for  the  X-Grid  (in  metres) 

1000. 

XEND  -  The  end  point  for  the  X-Grid  (in  metres) 

1000. 

XSTEP  -  The  step  size  in  X  (in  metres) 

0  . 

YSTART  -  The  starting  point  for  the  Y-Grid  (in  metres) 

0  . 

YEND  -  The  end  point  for  the  Y-Grid  (in  metres) 

2000  . 

YSTEP  -  The  step  size  in  Y  (in  metres) 

10  . 

ZSTART  -  The  starting  point  for  the  Z-Gnd  (in  metres! 

o . 

ZEND  -  The  end  point  for  the  Z-Grid  (in  metres! 

480  . 

ZSTEP  -  The  step  size  in  Z  (in  metres) 

10. 

DEPTH  -  The  depth  of  the  ocean 

480  . 

DENFIL  -  Name  of  the  file  containinq  the  density  profile  info. 

dawscn . pro 

NUMDEN  -  Number  of  points  in  the  density  profile,  DENFIL 

97 

DBGFIL  -  Name  of  the  file  to  write  the  debugging  information  to 

debug . dat 

VELFIL  -  Name  of  file  to  receive  output  velocity  component  data 

arrow.dat 

INTFIL  -  Name  of  the  file  to  receive  integration  data 

it.c  .  dat 

COMPU  -  Compute  the  U  component  velocity  ??  (Y)  or  (N) 

N 

COMPV  -  Compute  the  V  component  velocity  ??  (Y)  or  (N) 

Y 

COMPW  -  Compute  the  W  component  velocity  ??  (Y)  or  (N) 

Y 

Figure  2  :  Sample  Parameter  File 
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NUMTHETA  : 

MODE  : 

MAXTHETA : 
SUBDTH  : 

SUBSPD  : 

SUBLEN  : 
SUBRAD  : 
XSTART  : 

XEND  : 

XSTEP  : 

YSTART  : 

YEND  : 

YSTEP  : 


An  integer  in  the  range  2.. 30.  NUMTHETA  describes  the  number  of 
theta  values  to  divide  the  interval  0mo..7t/2  into.  The  amplitude 
functions  are  then  calculated  at  those  theta  values. 

An  integer  in  the  range  1..10.  MODE  specifies  the  internal  wave 
mode  that  we  wish  to  calculate  the  data  for. 

The  maximum  theta  value  allowed.  MAXTHETA  <  jt/2. 

A  real  in  the  range  0.. Depth  of  the  Ocean.  SUBDTH  specifies  the 
depth  of  the  submarine  in  metres.  The  Depth  of  the  Ocean  is  the 
input  parameter  DEPTH. 

A  real  in  the  range  O..Maximum  Speed.  SUBSPD  specifies  the  speed 
of  the  submarine  in  metres/second.  Maximum  speed  is  set  at  30  ra/s. 

A  real.  SUBSPD  specifies  the  length  of  the  submarine  in  metres. 

A  real.  SUBRAD  specifies  the  radius  of  the  submarine  in  metres. 

A  real.  XSTART  specifies  the  point  in  the  X-Plane  at  which  the 
integrations  will  be  started  from. 

A  real.  XEND  specifies  the  point  in  the  X-Plane  at  which  the 
integrations  will  end.  XEND  >  XSTART  >  0 

A  real.  XSTEP  specifies  the  size  of  the  X-Grid  (i.e.  the  spacing 
between  consecutive  points  in  the  X-Plane  at  which  an  integration 
will  be  calculated.  XSTEP  >  0 

A  real.  YSTART  specifies  the  point  in  the  Y-Plane  at  which  the 
integrations  will  be  started  from. 

A  real.  YEND  specifies  the  point  in  the  Y-Plane  at  which  the 
integrations  will  end.  YEND  >  YSTART  >  0 

A  real.  YSTEP  specifies  the  size  of  the  Y-Grid  (i.e.  the  spacing 
between  consecutive  points  in  the  Y-Plane  at  which  an  integration 
will  be  calculated.  YSTEP  >  0 
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ZSTART  :  A  real.  ZSTART  specifies  the  point  in  the  Z-Piane  at  which  the 

integrations  will  be  started  from,  (ie  first  depth  position) 

ZSTART  must  be  a  multiple  of  :  DEPTH  /  ((NUMDEN- 1  )/2) 

ZEND  :  A  real.  ZEND  specifies  the  point  in  the  Z-Plane  at  which  the 

integrations  will  end  (ie  last  depth  position). 

ZEND  must  be  a  multiple  of  :  DEPTH  /  ((NUMDEN- 1  )/2) 

ZSTEP  :  A  real.  ZSTEP  specifies  the  size  of  the  Z-Grid  (i.e.  the  spacing 

between  consecutive  points  in  the  Z-Plane  at  which  an  integration 
will  be  calculated. 

ZSTEP  must  be  a  multiple  of  :  DEPTH  1  ((NUMDEN- 1)/2) 

DEPTH  :  The  depth  of  the  ocean  as  described  by  the  density  profile.  DEPTH  >C 

DENFIL  :  A  character  string  containing  the  name  of  the  file  which  holds  the 

density  profile. 


**** 

**** 


DENSITY  PROFILE  MUST  HAVE  AN  ODD  **** 
NUMBER  OF  POINTS  **** 


NUMDEN  :  The  number  of  points  in  the  density  profile,  DENFIL. 

DBGFIL  :  A  character  spring  containing  the  name  of  the  file  to  receive  the 

debugging  information. 

VELFIL  :  A  character  string  containing  the  name  of  the  file  to  receive  the 

velocity  component  data. 

INTFIL  :  A  character  string  containing  the  name  of  the  file  to  receive  the 

integration  data. 

COMPU  :  A  character  indicating  whether  the  U  component  velocity  should  be 

calculated.  "Y"  or  "y"  indicates  "YES!,  Calculate  the  U  component" 
"N”  or  "n"  indicates  "No!,  Calculation  of  the  U  component  is  not 
required". 

COMPV  :  A  character  indicating  whether  the  V  component  velocity  should  be 

calculated.  "Y”  or  "y"  indicates  "YES!,  Calculate  the  V  component" 
"N"  or  "n"  indicates  "No!,  Calculation  of  the  V  component  is  not 
required". 

A  character  indicating  whether  the  W  component  velocity  should  be 
calculated.  "Y"  or  "y"  indicates  "YES!,  Calculate  the  W  component" 
"N"  or  "n”  indicates  "No!,  Calculation  of  the  W  component  is  not 
required". 


COMPW  : 
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2.2.  Density  Profile 


The  density  profile,  p(z)  is  contained  in  the  file  specified  in  die  parameter  DENTIL.  The  values  of 
p(z)  are  used  to  calculate  p(z)  according  to 


p  ( z ) 


p  '{z) 
p  ( z ) 


where  z  is  the  depji  ranging  from  0  (ie  the  ocean  surface)  to  DEPTH  (ie  the  ocean  floor) 
(NB  :  p'(z)  is  approximated  in  the  program  by  a  finite  difference  method) 


Figure  3  shows  a  sample  density  profile. 

NB  :  The  Density  Profile  must  have  the  following  format  : 

4  Header/Comment  lines  :  Gives  information  such  as  the  files  purpose  and  date. 

"NUMDEN"  Data  Lines  :  These  lines  contain  the  density  data.  TTie  first  line 

contains  the  density  at  the  Ocean  Floor  and  the 
NUMDEN1*’  Line  contains  the  density  of  the  Ocean 
Surface.  NUMDEN  has  to  be  odd  and  there  should 
be  no  blank  lines  at  the  end  of  the  file 


The  profile  shown  is  taken  from  T.W.  Dawson’s  DREP  Internal  Wave  Normal  Mode  -  Theoretical 
Background  Report,  1988.  [DREP  TM88-7] 


Dawson  Density  Premie, 
every  5  m. 

Depth  of  the  Ocean  =  480 

Number  of  Points  -  97 

1.0029176 

1.0029131 

1.0029087 

1.0029042 

1.0028996 

1.0028951 

1 .0028907 

1 .0028896 


Tony  Degg  7Z7J771 interpolated  ati 


1 .0000362 
1.0000062 
1.0000007 
1.0000000 

<EOF> _ 

Figure  3  :  Simple  Density  Profile 


NB  :  The  number  of  points  in  the  density  profile,  NUMDEN,  is  linked  implicitly  with  the 
depth  of  the  ocean,  DEPTH.  (NUMDEN-l)/2  should  divide  DEPTH  evenly 
i.e.  (  DEPTH  modulo  (NUMDEN-l)/2) )  =  0 


Internal  Waves  in  Submarine  Wakes 


8 


3.  Output 


The  program  produces  an  intermediate  output  file  containing  the  integration  data,  but  the  main 
output  file  contains  the  velocity  component  data. 

3.1.  Velocity  Component  Output  File 


The  velocity  component  output  file  contains  a  stream  of  six  (6)  columns  of  numbers  which  can  be 
used  by  MATLAB  (or  other  graphing  utilities)  to  create  several  different  type  of  graphs.  Figure  4 
shows  a  sample  velocity  component  output  file  : 


100  . 

0. 

0  . 

0. 

0  .  OOOOOE-OO 

-6 . 92253E-07 

100. 

10. 

0  . 

0. 

7 . 92248E-04 

-5 . 79579E-07 

100. 

20. 

0. 

0. 

1 . 17258E-03 

-3  . 123 13E-07 

100  . 

30  . 

0  . 

0  . 

1 . 1 1829E-Q3 

-4 . 18162E-08 

100. 

40  . 

0. 

0. 

8 . 22264E-04 

1 .45645E-07 

100. 

50. 

0  - 

0. 

4 . 7 3  894E-04 

2.32705E-07 

100  . 

SO  . 

0. 

0  . 

1 . 83  084E-04 

2 . 40433E-07 

100  . 

70. 

0 « 

0. 

- 1 . 63455E-05 

2.07768E-07 

100. 

so. 

0  . 

0 

-1 .41567E-04 

1.55930E-07 

100. 

90. 

0 . 

0 

-2 . 05422E-04 

7 . 54386E-08 

100. 

100. 

0. 

0. 

-2 . 11459E-04 

8 . 57241E-08 

200  . 

0. 

0 . 

0. 

0. OOOOOE-OO 

-6 . 96273E-07 

200. 

10. 

0. 

0. 

8 . 92248E-04 

-3 . 76574E-07 

200. 

20. 

0  . 

0  . 

1 . 47278E-03 

-2 . 143 16E-07 

200  . 

30. 

0. 

0 

1.14799E-03 

-4 . 58662E-08 

200. 

40. 

0. 

0 

8 . 52764E-04 

2 . 46645E-07 

200. 

50. 

0  . 

0. 

4 . 70934E-04 

4 . 63655E-07 

1000. 

SO. 

480  . 

0. 

1 . 41223E-0S 

-6 . 67  499E-07 

1000. 

70. 

0  . 

1 . 42418E-06 

-8 . 10197E-07 

Figure  4  :  Sample  Velocity  Component  Output  File 


It  consists  of  6  columns  of  data.  The  last  three  columns  correspond  to  the  integrals  of  U,V  and  W, 
respectively,  at  the  point  X,Y X  specified  in  the  first  three  columns.  (NB:  in  this  particular  run  of 
WAKE  only  the  V  and  W  velocity  components  were  calculated,  the  U  column  therefore  contains 
only  zeroes).  It  is  up  to  MATLAB  or  any  other  graphing  utility  to  separate  the  data  correctly  and 
produce  meaningful  graphs.  See  Section  4  :  MATLAB. 

Several  different  types  of  graphs  can  be  produced.  Figure  5  shows  a  2-D  Wave  Plot  and  a  3-D 
mesh  plot  of  the  V  Velocity  Component  taken  at  the  ocean  surface  with  the  X-Y  grid  extending  6 
and  2  kilometres  respectively.  Figure  6  shows  a  Matrix  Intensity  Plot,  produced  using  PLOTZ,  XV 
and  MATLAB.  Figure  7  shows  the  V  and  W  Component  Velocity  Contour  Plot.  Figure  8  shows  a 
zoomed  in  image  of  Figure  7  with  arrows  indicating  the  field  flow. 
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X  Distance  in  Metres 


Figure  5  :  3-D  Mesh  and  2-D  Wave  Plots 
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Figure  6  :  Matrix  Intensity  Plot 


Depth,  in  metres  Depth,  in  metres 
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3.2.  Integration  File 


The  intermediate  output  produced  by  WAKE  is  the  integration  file.  This  file  contains  the  data 
required  for  the  integration.  It  consists  of  two  parts,  a  Header  and  a  Main  Body. 

Header  :  The  header  contains  some  of  the  parameters  used  in  calculating  the 

amplitude  functions.  These  are  : 

Mode  Number,  Number  of  Depth  Values,  Depth  Step, 

Number  of  0  Values,  Submarine  Depth  and  Speed. 

Main  Body  :  The  main  body  contains  for  each  0,k  pair  the  amplitude  functions 

A(0)  and  B(0)  for  all  depth  values. 


Figure  9  shows  a  sample  integration  data  file  : 


Mode  Number 

Number  c i  Z  Values 

49 

Z  St ec 

.ccoooooooooc: 

Number  cf  Theca  Values 

30 

Theca  Seep 

:.:S0138229491$:-02 
Submarine  Soeed 

2.5000000000000 
Submarine  Depth 

so.ooooooooooo: 

HEADS?. 

Contains  me  parameters 

used  to  generate  me 

Am.pl it'-de  Functions 

Theca 

C  -95S7  54-;  5794732 

Wave  Number,  K 

4 .99032120121353-13 

A (Theta) 

-2 .54190911781052-20 
-2.54137792918353-20 
-2.53977661974263-20 
-2  .53709788888993-20 

B(Theta) 

-2 . 5419091 17310SD- 20 

2 . 1288060573866D-G7 

4  .29660423382443-07 

6. 4307667182042 D- 07 

2 

-480.00000000000 

-470.00000000000 

-460.00000000000 

-450.00000000000 

MAIN  BODY 

Repeated  NUM THETA  times 

. 

from  $„>n. .  x/2  -8,^; 

-3  .22281720031933-21 
-1 . 59487475686453-21 
3.45973986536623-23 

3.24970131603363-05 

3  . 2 650882233180D-05 

3. 2652830218816 D- 05 

-20.000000000000 

-lO.OOOCCCOOOCCO 

O.OOOOOOOCOCCO 

Theta 

1.5502949445000 

Wave  Number ,  K 

0.65904859979995 

A (Theta) 

1  . 984295891 144S-10I 

3 .05384485606433-99 
6.049225025  62  96 D- 97 

B(Theta) 

1.9842958911445-101 

5 . 5573161700118 D- 99 

1 .10519148475123-96 

2 

-480.0000000000C 

-470.00000000000 

-460.00000000000 

-2.76223226063833-07 
-2 . 5742411 59 5178D-09 
-1.31204182739293-11 
4.27170128745243-17 

2  .43459127650423-07 

4. 07509629277423-09 

2 .43720369709 OSD- 11 

2 .41837726013753-13 

-30.000000000000 

-20.000000000000 

-10.000000000000 

0.000000000000 

Figure  9  :  Sample  Integration  Data  File 
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4,  MATLAB 

To  generate  the  graphs  shown  in  Figures  5, 6,7,8  requires  a  graphics  package.  The  one  described 
here  is  called  MATLAB  and  can  be  run  both  on  PC’s  and  on  SUN  Workstations. 

To  get  started,  invoke  the  matlab  package  with  the  command  : 

>  matlab 

After  a  moment,  the  MATLAB  banner  and  a  prompt  like  "»"  will  appear.  The  MATLAB 
interpreter  is  awaiting  instructions.  Load  in  the  velocity  component  data  file  generated  by  WAKE 
as  follows  : 

» load  <filename> 

After  the  file  has  been  loaded,  End  out  the  size  of  the  matrix  that  MATLAB  has  loaded  the  data 
into  by  executing  the  command  : 

»  whos 

MATLAB  will  return  with  an  answer  that  looks  something  like  : 
doiphin>  matlab 

<PRO-MATLAB> 

(c)  Copyright  The  MathWorks,  Inc.  1984-1991 
All  Rights  Reserved 
Version  3.5i  18-Jul-1991 

HELP,  DEMO,  INFO,  and  TERMINAL  are  available 

>> 

>>  load  output.dat 
>>  whos 

Name  Size  Total  Complex 

output  5015  by  6  30090  No 

Grand  total  is  (30090  *  8)  =  240720  bytes, 

>> 

Figure  10  :  Finding  out  about  current  variables  and  their  sizes 
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Notice  the  matrix  is  a  six  column  matrix.  To  generate  diagrams  like  Figure  5  and  Figure  6  requires 
reshaping  the  matrix  into  an  appropriate  format.  This  is  done  as  follows  : 

»  u  =  reshape(filename(:,4),m,n)  ; 

»  v  =  reshape(fllename(:,5 ),m,n)  ; 

»  w  =  reshape(filename(:,6),m,n)  ; 


where  : 

FILENAME  =  Name  of  the  data  file  generated  by  WAKE  without  its  file  extension 
The  symbol  means  "all".  Thus.  A(:,4)  means  extract  all  rows  of 
the  fourth  column  of  the  matrix  A,  whereas  A(4,:)  means  extract  all 
columns  of  the  fourth  row  of  the  matrix  A. 

M  =  Number  of  Y-Grid  points  generated  by  WAKE 

(i.e.  M  =  ((YEND- Y START)/Y STEP)  +  1 
N  =  A  number.  N  is  the  number  of  X-Grid  points  generated  by  WAKE 

(i.e.  N  =  ((XEND-XSTARTVXSTEP)  +  1 
U  =  U  is  the  resultant  M  x  N  matrix  corresponding  to  the  U 

component  velocity  data  (if  COMPU  =  'Y'  ie  WAKE  has  generated 
the  data  for  the  U  velocity  component) 

V  =  V  is  the  resultant  M  x  N  matrix  corresponding  to  the  V 

component  velocity  data  (if  COMPV  =  ’Y’  ie  WAKE  has  generated 
the  data  for  the  V  velocity  component) 

W  =  W  is  the  resultant  M  x  N  matrix  corresponding  to  the  W 

component  velocity  data  (if  COMPW  =  'Y'  ie  WAKE  has  generated 
the  data  for  the  W  velocity  component) 

4.1.  3-D  Mesh  Plots 


To  generate  a  3-D  Mesh  Plot  of  the  V  Velocity  Component  as  shown  in  Figure  5.  execute  the 
following  command  : 

»  mesh(v,[hr,e]) 


where  : 

V  =  MxN  matrix  containing  the  V  velocity  component  data 

HR  =  Horizontal  Rotation 

E  =  Elevation 

NB  :  A  viewpoint  matrix  [HR.E]  of  [80,80]  is  recommended  as  a  start. 

A  title  may  be  put  on  the  plot  as  follows  : 


» title(’3-D  Mesh  Plot’) 
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4.2.  2-D  Wave  Plots 


To  generate  a  2-D  Wave  Plot  of  the  V  Velocity  Component  as  shown  in  Figure  5.  execute  the 
following  commands  : 

First  generate  an  offset  matrix  as  follows  : 

»  offset  =  ones(M,l )*Unspace(0,  TOP, N)/l 00000.  ; 

where 

N  =  Number  of  columns  in  our  matrix  V 

M  =  Number  of  rows  in  our  matrix  V 

TOP  =  Arbitrary  scale  for  the  Y-Axis  of  the  plot. 

Next  generate  a  scaling  matrix  as  follows  : 

»  scale  =  linspace(ymin,ymax,M )  ; 


where 

N  =  Number  of  columns  in  out  matrix  V 

YMIN  =  Starting  value  of  the  Y-Axis 

YMAX  =  End  value  of  the  Y-Axis 

Now  define  a  new  matrix,  X  which  is  the  matrix  V  with  each  column  offset  slightly  as  follows  : 

»  x  =  v  +  offset ; 

The  2-D  Wave  Plot,  with  an  appropriate  title,  is  generated  in  the  Graphics  Window  with  the 
following  commands  : 

»  plct(x, scale,  ’b-’) 

» title(’2-D  Wave  Plot’) 

4.3.  Matrix  Intensity  Plots 

To  generate  a  Matrix  Intensity  Plot  of  the  V  Velocity  Component  like  Figure  6.  requires  a  package 
called  PLOTZ.  First  save  the  matrix  A  to  a  file.  This  is  done  by  executing  the  command 

»  save  wave  v 

This  generates  a  file  called  WAVE. MAT,  in  the  current  directory,  which  contains  the  data  from  the 
matrix  A.  Now  exit  from  MATLAB  with  the  command 


»  quit 
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To  start  up  PLOTZ  and  generate  the  matrix  intensity  plot  execute  the  following  command: 
>  plotz  wave. mat  [1,2, 3, 4]  xmin  xmax  ymin  ymax 


where 


WAVE.MAT  =  File  generated  by  MATLAB  containing  the  velocity  component  data 
[1, 2,3,4]  =  Choose  one  of  the  following  options  : 

1  =  Grey 

2  =  Colour 

3  =  4  Tones  of  Green 

4  =  Black  &  White  (default) 

XMIN  =  Starting  value  of  the  X-Axis 

XMAX  =  End  value  of  the  X-Axis 
YMIN  =  Starting  value  of  the  Y-AXIS 
YMAX  =  End  value  of  the  Y-AXIS 


4.4.  Contour  and  Field  Flow  Plots 


To  generate  the  contour  and  quiver  plots  as  shown  in  Figure  7  and  Figure  8,  run  WAKE  fixing  X 
at  some  point  behind  the  submarine,  generating  the  V  and  W  component  velocities  in  the  Y-Z 
Plane.  First  load  the  data  and  reshape  the  matrix  into  an  appropriate  format  as  follows  : 

» load  <filename> 

»  v  =  reshape(filenam>  • ,5),m,n )  ; 

»  w  -  reshape(filename(:,6),m,n)  ; 

where  : 


FILENAME  =  Name  of  the  data  file  generated  by  WAKE  without  its  file  extension 
M  =  Number  of  Y-Griu  points  generated  by  WAKE 

(i.e.  M  =  ((YEND-YSTARTVYSTEP)  +  1 
N  =  A  number,  N  is  the  number  of  Z-Grid  points  generated  by  WAKE 

(i.e.  N  =  ((ZEND-ZST ART)/ZSTEP)  +  1 

V  =  V  is  the  resultant  M  x  N  matrix  corresponding  to  the  V  component 

velocity  data  (if  COMPV  =  ’Y’  ie  WAKE  has  generated  the  data  for 
the  V  velocity  component) 

W  =  W  is  the  resultant  M  x  N  matrix  corresponding  to  the  W  component 

velocity  data  (if  COMPW  =  ’Y’  ie  WAKE  has  generated  the  data  for 
the  W  velocity  component) 
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4.4.1.  Contour  Plots 


Set  the  scale  on  the  Y  and  Z  axis  to  correspond  to  the  data  generated. 

»  yscale  =  linspace(ystart,yend,M)  ; 

»  zscale  =  linspace(-zend,zstart,N); 

Plot  and  keep  the  empty  axis  as  follows  : 

»  plot([y start  yend],  [-zend  zstartj,’.’) 

»  hold  on 

Plot  the  V  component  velocity  contour  as  follows  : 

»  vt  =  v’; 

»  contour(vt,y scale, zscale,  ’b-’) 

VT  is  now  the  transpose  of  the  matrix  v.  This  is  required  to  get  the  correct  orientation. 
Appropriate  titles  and  axis  labels  are  added  with  the  following  commands  : 

»  titlef’V  Component  Velocity  Contour  Plot,  at  1000m  behind  the  Submarine’) 

»  xlabel  (’ Distance  Y,  in  metres) 

»  ylabel  (’Depth,  in  metres) 

4.4.2.  Field  Flow  Plots 


To  produce  a  contour  plot  with  field  flow  arrows  overlayed,  display  in  the  Graphics  Window  the 
contour  plot  and  execute  the  following  commands  : 

»  hold  on 
»  quiver(vt,wt) 

where  VT  and  WT  are  the  transposes  of  the  V  and  W  velocity  component  matrices.  This  will 
produce  a  contour  plot  with  field  flow  arrows  overlayed.  The  arrows  are  scaled  so  that  they 
indicate  not  only  direction  but  magnitude.  The  scaling  of  the  arrows  become  too  small  if  we  are 
looking  at  a  wide  area,  hence  it  is  necessary  to  look  at  small  portions  of  the  plot  at  any  one  time  to 
see  the  structure  of  the  field  flow. 
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To  look  at  a  section  of  the  plot  like  Figure  8  execute  the  following  : 

»  clg 
»  hold  off 

»  contour(vt(zindstart:zindend,yindstart:yindend), ’b-’) 
»  hold 

»  quiveri vt(zindstart:zindend,yindstart:yindend), 

wt(zindstart:zindend,yindstart:yindend)) 
»  title(’V  Contour  Plot  with  Field  Flow) 


Where  : 

ZINDSTART,ZINDEND,YINDSTART,YINDEND  define  the  submatrix  of  the 
matrix  VT. 

This  has  the  effect  of  zooming  in  on  a  portion  of  the  plot  so  that  we  can  see  the  structure  of  the 
field  flow. 

4,5.  Printing 

To  print  within  MATLAB.  display  the  graph  required  in  the  graphics  window  and  then  execute  the 
following  commands  : 

»  meta  <fdename> 

This  creates  a  file  called  FILENAME.MET  in  the  current  directory.  Create  a  postscript  file  of  the 
plot  as  follows  : 

»  Igpp  filename.met  -dps  -ol 

This  generates  a  postscript  file  of  the  plot  which  can  be  sent  to  a  postscript  printer  in  the  usual 
manner. 

For  more  information  on  MATLAB  and  MATLAB  commands  consult  the  MATLAB  user  manual. 
To  print  a  matrix  intensity  plot  like  Figure  6,  save  the  plot  generated  by  PLOTZ,  as  follows  : 

1) .  Click  on  the  plot  with  the  right  button  of  the  mouse 

2) .  Select  Save  to  Raster 

3) .  Quit  out  of  PLOTZ  by  clicking  again  on  the  plot  with  the  right  mouse  button  and 

selecting  Quit 
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PLOTZ  has  now  generated  a  Raster  file  called  "PLOTZ.RAS".  To  print  this  file  convert  it  to 
PostScript  format  as  follows  : 

1) .  Execute  the  Raster  file  viewing  package  XV  as  follows  : 

>  xv  plotz.ras 

2) .  Click  with  the  right  button  of  the  mouse  on  the  plot.  This  brings  up  a  menu  with 

various  options  and  parameters. 

3) .  Select  Save.  This  brings  up  another  menu  with  more  options. 

4) .  Select  PostScript  and  rename  the  output  file  (if  desired)  and  click  on  OK 

5) .  This  brings  up  a  final  menu  with  all  the  PostScript  options  (i.e.  portrait,  landscape, 

paper  size,  etc).  Select  the  parameter  desired  and  click  on  OK 

6) .  A  rotating  fish  indicates  that  XV  is  converting  the  Raster  file  to  a  Postscript  file. 

7) .  Quit  out  of  XV  and  send  the  PostScript  file  just  generated  to  a  PostScript  printer  in 

the  usual  manner. 
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To  compile  the  WAKE  program  the  following  files  need  to  be  copied  into  a  working  directory. 

Access  to  the  IMSL  library  Routines  DQDAG,  CSAKM,  CSVAL  (fortran  version)  is  also  required. 

MAKEFILE  :  Compiles  all  the  routines  into  a  library  and  links  them  all  together. 

Figure  12  shows  a  sample  makefile. 

COMP  :  A  script  file  which  compiles  the  source  fortran  files  to  object  files  and 

creates/updates  a  library  called  WAKE  which  contains  all  the 
necessary  routines  and  functions.  Figure  1 1  shows  a  sample  compile 
file. 

DENSITY.DAT  Name  of  the  file  which  contains  the  density  profile. 

WAKE.PAR  Contains  all  of  the  parameters  WAKE  requires 

YVAKE.F  Main  Program  controls  all  the  following  subroutines 

BOUND.F  Searches  for  changes  of  sign  in  the  function  D(9Jt)  for  a  given  0  and 

saves  the  intervals  in  a  matrix 

COMPUTE_A.F  Computes  the  Amplitude  Functions  A(8),  B(9). 

COMPUTE_H.F  Converts  a  real  depth  value  into  an  integer  grid  point. 

DENSITY.F  Reads  in  the  density  profile  and  computes  MU(z). 

DET.F  Returns  the  value  of  D(o,k),  where  D  =  Determinant  of  the 

Wronskian  of  the  ODE. 

EXPANALYT.F  Calculates  the  analytical  solutions  for  an  exponential  density  profile. 

FINDINTERVAL.F  Finds  an  interval  within  which  a  root  of  the  function  passed  occurs. 

FIND_ZEROES.F  Returns  the  0(0)  values. 

GETSIGVALS.F  Uses  a  quadratic  mapping  technique  to  obtain  an  array  of  o 
values  biased  at  the  endpoints. 

INIT.F  Reads  in  the  initial  parameters  from  the  parameter  file. 

INTEGRA TE.F  Performs  the  integrations  using  the  IMSL  routine  DQDAG. 

LINEAR.F  Performs  a  linear  extrapolation  to  obtain  a  guess  for  the  value  of  K  at 

the  current  o  value  using  two  previous  (0,k)  pairs. 

LOOKUP.F  Looks  up  the  interval  matrix  to  obtain  the  intervals  where  roots  of  the 

function  D(0.k)  occur. 

READINT.F  Reads  in  the  integration  data  for  a  fixed  depth. 

RTBIS.F  Numerical  Recipes  Bisection  Method. 

RTSEC.F  Numerical  Recipes  Secant  Method. 

UINT.F  Computes  the  U  velocity  integrand. 

VINT.F  Computes  the  V  velocity  integrand. 

WINT.F  Computes  the  W  velocity  integrand. 

WRITEHEAD.F  Writes  the  header  information  to  the  integration  file. 

WRITEINT.F  Writes  the  integration  data  to  the  integration  data  file. 

YVFUNCNS.F  Computes  the  functions  W,,  W/,  W,,  W/. 


COMMON.BLK  Contains  all  of  WAKE’s  common  variables. 


U)  XI 
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6.  Program  Description 


The  program  can  be  divided  into  four  main  stages  : 


STAGE  1 
STAGE  2 
STAGE  3 
STAGE  4 


Initialisation 

Binding  the  roots  of  the  function  D(a,k) 

Solve  the  O.D.E.  and  Calculate  the  Amplitude  Functions  A(0),  B(0) 
Form  the  3  direction  (x,y,z)  integrands  U,V,W  and  integrate  them 


These  stages  are  described  in  the  following  sections. 
6.1.  Staee  1  :  Initialisation 


Stage  1  sets  up  the  initial  conditions.  Figure  13  descnbes  diagrammahcally  the  structure  for 
Stage  1 


The  initial  parameters  are  first  read  in  from  the  parameter  file.  WAKE.PAR,  see  Section  2.1  for  a 
description  on  the  input  parameter  file.  Next  the  program  reads  in  from  the  file  specified  in  the 
parameter,  DENFIL,  the  density  profile,  see  Section  2.2  for  a  description  of  the  density  profile. 


6.2.  Stage  2  :  Binding  the  roots  of  DlQ.k-) 

Once  the  initial  parameters  and  density  profile  has  been  read  in.  consider  the  dispersion  relation 
fsee  E.O.  Tuck  :  Appendix  3  Submarine  Internal  Waves.  1992]. 

"  Our  concern  here  is  very  much  with  dispersion  relations  for  internal  waves,  namely  relations 
between  wave  speed  c  and  wave  number  k.  For  convenience,  we  use  instead  of  c  a  quantity 
proportional  to  its  reciprocal  square,  namely 


where  g  is  gravity  Then  we  need  a  connection  between  k  and  a,  e.g.  k  =  K(o).  One 
of  our  first  tasks  is  to  compute  this  relation  for  a  given  density  distribution."  [EOT  92] 
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Figure  14  :  Dispersion  Relations 


Figure  14  shows  the  dispersion  relation  for  the  first  4  internal  waves  (The  unseen  line  6=k  would 
correspond  to  the  surface  wave).  We  can  see  that  as  9  —>  n/2  the  graphs  get  very  close  together 
and  at  0  =  rc/2  there  are  in  fact  an  infinite  number  of  internal  waves  (this  becomes  a  problem  later 
when  the  program  tries  to  pick  out  one  particular  internal  wave  and  stay  within  that  one  mode). 
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CALCULATE  0  POINTS 

USING  QUADRATIC  MAPPING 

MHMM 

\ 

! 

Figure  15  :  Stage  4  :  Binding  the  roots  of  D(0,k) 


The  dispersion  relation  needs  to  be  accurately  determined  for  multiple  modes.  The  first  step  is  to 
accurately  determine  the  first  mode  as  this  acts  as  an  upper  bound  on  the  other  modes.  The  first 
mode  is  determined  in  the  following  manner  : 

Each  mode  has  a  9min,  denoted  0'node  (e.g.  the  0min  value  for  the  first  mode  is  denoted  01).  Below 
this  0min  value  internal  waves  do  not  exist  for  that  mode  (therefore  below  the  0min  value  for  the  first 
mode,  no  internal  waves  exist).  The  interval  0‘..THETAMAX  is  divided  into  NUMTHETA 
subintervals,  each  of  size  0SIep.  The  values  of  K  for  0‘  and  0'+0s(ep  are  determine  by  an  exhaustive 
search. 
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Once  these  first  two  values  have  been  determined  a  linear  extrapolation  guess,  K^,  for  the  value 
of  K  at  the  next  0  value  can  be  obtained. 

[  It  was  found  that  perturbations  in  the  dispersion  curve  resulted  in  the  linear  extrapolation 
projected  below  the  first  mode  curve.  To  overcome  this  a  test  on  the  sign  of  the  function 
D(0,k)  at  the  point  is  performed.  If  the  function  is  negative  then  the  extrapolation 

must  be  above  the  first  mode,  as  desired  (the  probability  of  the  extrapolation  projecting  into 
the  third  mode  region  where  the  function  is  also  negative  is  considered  to  be  negligible).  If 
the  function  is  positive  then  the  extrapolation  is  below  the  first  mode  curve  so  the  previous 
(0Jc)  pair  is  used  for  the  extrapolation.  This  process  of  using  progressively  previous  (0,k) 
pairs  is  repeated  until  a  satisfactory  linear  extrapolated  guess,  has  been  obtained  ] 

The  process  is  continued  from  0‘+0slep  through  to  THETAMAX,  at  each  step  using  the  two 
previous  (0,k)  pairs  to  obtain  a  linear  extrapolated  guess,  Kg^,  for  the  value  of  K  at  the  next  0 
value.  Once  has  been  obtained,  the  program  searches  downwards  in  steps  of  looking  for 
intervals  where  a  sign  change  in  the  function  D(0,k)  occurs  (i.e.  searches  for  roots  of  the  function 
D(0,k)  )  For  any  given  0  value,  the  minimum  number  of  sign  changes  that  should  occur  is  known 
through  the  determination  of  the  0mode  values.  If  the  number  of  sign  changes  found  is  less  than  the 
number  of  sign  changes  that  should  occur  then  the  K^,,  value  used  must  have  been  too  large,  so 
K^p  is  halved  and  the  search  for  sign  changes  is  repeated  until  the  correct  number  of  sign  changes 
that  should  occur  have  been  found.  Once  the  correct  number  of  sign  changes  have  been  found  the 
intervals  are  saved  in  a  matrix  with  the  columns  corresponding  to  the  various  0  values  and  the  rows 
containing  the  intervals  where  sign  changes  occur. 

When  the  program  wants  to  know  the  K  value  for  a  certain  mode  at  a  given  0  value  all  that  is 
required  is  to  look  up  in  the  INTERVAL  matrix  to  get  the  interval  where  the  sign  change  occurs 
for  the  given  mode  and  use  the  Bisection  Method  to  hone  in  on  the  root.  If  a  higher  mode  than  the 
first  mode  is  selected  then  the  program  still  calculates  NUMTHETA  0  values  between  0modc  and 
THETAMAX,  however,  10  preliminary  calculation  are  performed  from  01  to  0mode. 
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6.3.  Stage  3  :  Solving  the  ODE 

Figure  16  shows  diagrammatically  the  structure  for  stage  3 


Figure  16  :  Stage  3  Solving  the  ODE 


In  this  stage  we  solve  the  ODE 

(p&’V-  (i^p  +  opOw  =  0 

and  calculate  the  amplitude  functions  A(0),  B(0). 

First  we  divide  the  interval  9min..7i/2  into  NUMTHETA  intervals  and  calculate  0slep, 

*  A 
~Z  Vain 

0  —  ^ 
step  ~  NUMTHETA 


The  Amplitude  functions  are  then  calculated  at  0  =  0min  ..  (Jt/2-0sttp)  in  steps  of  0slcp. 

So  for  each  0  we  calculate  the  corresponding  a  value. 

o  =  Ksec2 (9) 
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By  fixing  c  and  varying  k  we  now  find  an  interval  k{..k^  where  D(a,k)=0  with  k,  <  b  <  k,.  Once 
this  interval  has  been  found  (  i.e.  the  root  has  been  bounded)  we  use  the  Secant  Method  to  hone  in 
on  the  root.  At  this  stage  we  have  a  0(a)  and  corresponding  k  pair,  we  use  this  pair  to  calculate 
the  Amplitude  functions  A(0),  B(0)  as  shown  in  Section  6.3.1. 


6.3.1.  Amplitude  Functions  A(0).  B(0) 


For  a  detailed  description  of  the  method  of  calculating  A(0)  and  B(0)  see  [E.O.  Tuck  1992  : 
Submarine  Internal  Waves  :  Appendix  1] 

NB:  The  amplitude  functions  are  not  defined  explicitly  in  EOT  92,  they  appear  implicitly  in 
equation  (4.1) 


A  (0) 


kw0 

dD(Q,k) 

dk 


B(Q ) 


dP(Q,k) 

dk 


where 


|  -tfi(h)  Wx(z)  :  Zih 
|-ft^(h)W2U)  :  hzzsQ 


and  h  =  Submarine  Depth 

W  =  W,(z)  and  W  =  W,(z)  are  two  separate  solutions  of  the  ODE 

<pWV  -  <*2p  +  op')  W  =  0 


where  a  =  xsec2©  ;  p  =  Density 

Once  the  Amplitude  Functions  have  been  calculated  the  data  is  written  to  the  data  file  specified  in 
the  parameter  INTFIL  as  specified  in  the  parameter  file  "WAKE.PAR". 


Internal  Waves  in  Submarine  Wakes 


29 


6.4.  Stage  4  :  Inteerands  and  Integration 

In  this  stage  of  the  program  we  form  the  directional  velocity  component  integrands  U,  V,  Z  and 
integrate  them.  Figure  17  shows  diagrammatically  the  structure  for  Stage  4 


The  integration  is  performed  by  the  IMSL  package  routine  DQDAG.  This  routine  expects  a 
continuous  function  (we  achieve  this  by  taking  our  discrete  vectors  which  describe  the  functions 
A(0),  B(0)  and  K(0)  and  performing  Akima  Cubic  Spline  Interpolations,  using  IMSL  routines 
DCSAKM  and  DCSVAL  as  required  to  obtain  continuous  smooth  functions)  and  an  interval  to 
integrate  over  (in  our  case  0mm..rc/2).  The  IMSL  integration  routine  then  requires  some  desired 
accuracy  and  quadrature  rules  to  be  set.  DQDAG  is  a  general  purpose  integrator  that  uses  a 
globally  adaptive  scheme  in  order  to  reduce  the  absolute  error.  It  subdivides  the  integrating 
interval  and  uses  a  (2k+l)-point  Gauss-Kronrod  rule  to  estimate  each  subinterval.  It  was  found  that 
a  Gauss-Kronrod  Rule  with  20-41  points,  and  a  requested  accuracy  of  l.D-7  produced  the  best 
results. 

The  amplitude  data  A(0),  B(0)  is  stored  in  the  data  file  INTFEL.  The  amplitude  data  is  read  back 
into  the  program  so  that  we  have  the  arrays  A,  B  containing  the  amplitude  functions  for  all  0 
values  at  a  specific  depth.  (The  file  INTFIL  contains  the  amplitude  functions  for  all  depths  at  a 
specific  0  value). 
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The  first  step  is  to  define  the  Z-Grid  over  which  the  integrations  will  take  place,  then  for  each 
depth  position  fie  each  point  in  the  Z-orid)  the  A(0),  B(9)  and  K(9)  is  read  in.  The  next  step  is  to 
define  the  X-Y  Grid  over  which  the  integrations  will  take  place.  For  each  point  in  the  X-Y-Z 
subspace  the  integration  is  calculated  for  each  of  the  requested  components  and  the  results  are 
written  to  the  output  velocity  data  file  as  described  by  the  parameter  VELFTL. 

The  U  velocity  component  corresponds  to  the  X  direction  (ie  in  the  same  direction  that  the 
submarine  is  moving).  The  V  velocity  corresponds  to  Y  direction  (ie  sideways  to  the  direction  that 
the  submarine  is  moving).  The  W  velocity  corresponds  to  the  Z  direction  (ie  the  depth  direction). 
The  following  three  equations  describe  the  integrals  of  the  three  velocity  components. 


w 

U 


2V 

it 


x/2 

J"  cos  (.fcxcosQ)  cos(kysinQ) 


'  sin(jrcos0L/2) 
L/2 


A  (8)  d6 


u 

U 


2V 

it 


x/2 

j*  sin(kxcos0)  cos(kysinQ) 

Alin 


‘  sin  (itc os8  L/2) 
L/2 


B{0)  cos0  dd 


v 

U 


2V 

it 


X/2 

j*  cos(Jcxcos0)  sin(icysin0) 


'  sin(Jccos8 L/2) 
L/2 


B(0)  3in0d0 


where  U 


Submarine  Speed 
Submarine  Volume 
Submarine  Length 

Position  in  the  X-Plane  relative  to  the  submarine 
Position  in  the  Y-Plane  relative  to  the  submarine 
Amplitude  Function,  A 
Amplitude  Function.  B 
Wave  Number 
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7.  Testing 


Stage  1  of  the  program  consists  of  reading  in  parameters  and  the  density  profile.  This  section  can 
be  tested  by  using  the  DBXTOOL  debugging  utility  and  stepping  through  each  line  of  the  code 
checking  that  the  correct  values  are  read  in  to  the  right  variables. 

Stage  2  of  the  program  is  concerned  with  binding  the  roots  of  the  function  D(QJc).  The  method 
finally  chosen  (i.e.  Linear  extrapolated  guess  and  then  downwards  search  for  changes  in  sign)  was 
only  the  last  of  many  techniques  investigated,  (others  included  searches  using  the  secant  method, 
fixing  either  K  or  6).  This  final  technique  proved  to  be  the  most  robust  especially  in  terms  of 
distinguishing  the  higher  order  modes.  This  section  of  the  code  was  tested  using  the  DBXTOOL 
debugging  utility  to  ensure  that  the  K  value  the  program  hones  in  on  does  indeed  lie  between  the 
interval  for  that  mode. 

Stage  3  of  the  program  is  concerned  with  solving  the  ODE  and  computing  the  Amplitude  Functions 
A(0)  and  B(0).  For  an  exponential  density  stratified  ocean,  there  exists  exact  analytical  formulae 
for  the  functions  W,.  W/.W,,  W;'.  the  solutions  of  the  ODE  and  the  wronskian.  D  as  follows  : 

Let 


X  =  Jdisc 


where 

disc  =  45o  -  8:  -  4kz 

5  =  ln(Density  Ratio),  the  natural  logarithm  of  the  ratio  between 

the  density  of  the  ocean  floor  to  the  ocean  surface 
k  =  Wave  Number 

a  =  Wave  Speed  Parameter 


where 
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therefore  our  set  of  equations  become: 

[( ll£lj 
K  -  |w,  ♦  ~2 —  \(2k-b)costiUl'z) )  -  xsL-M' k 

K^M^)  -  Hir)l 

iw‘t±r  -  w^-r)] 


These  formulae  are  contained  in  the  routine  "EXPANALYTIC".  If  the  program  is  run  with  an 
exponential  density  profile  the  program  produces  tables  which  can  be  used  to  verify  that  the 
analytic  and  numerical  solutions  agree.  Figure  18  and  Figure  19  are  tables  of  data  taken  from  the 
debugging  file  which  show  the  numerical  and  analytic  solutions  for  WQ,  W,,  W/,  W:,  W2;  and  the 
amplitude  functions  A(0)  and  B(8).  These  tables  indicate  that  the  numerical  and  analytic  solutions 
agree  to  a  satisfactory  number  of  significant  figures. 
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Stage  4  of  the  program  forms  and  integrates  the  three  velocity  components  U,V  and  W.  To 
integrate  these  velocity  components  using  the  IMSL  routine  DQDAG  requires  continuous  forms  for 
the  functions  A(0),  B(0)  and  K(0).  However  A(0),  B(0)  and  K(0)  exist  only  as  discrete  tabulated 
functions,  therefore  some  form  of  interpolation  was  required. 

Three  different  interpolation  techniques  were  investigated.  These  were  : 

1) .  Linear  Interpolation 

2) .  Natural  Cubic  Spline  Interpolation 

3) .  Akima  Cubic  Spline  Interpolation 

Figure  20  shows  the  three  different  interpolation  methods,  used  to  produce  a  continuous  function 
for  K(0). 

The  linear  interpolation  is  quite  good  and  would  be  satisfactory  for  K(0),  however  for  A(0)  and 
B(0),  it  was  suggested  that  a  smoother  interpolation  technique  would  be  required. 

The  next  technique  investigated  was  Natural  Cubic  Spline  Interpolation.  This  technique  had  some 
unforseen  side-effects.  One  of  the  properties  of  natural  piece-wise  cubic  splines  is  that  the  first 
derivative  at  the  end  of  one  interval  must  equal  the  first  derivative  at  the  start  of  the  next  interval, 
however  no  attempt  is  made  to  preserve  the  function’s  shape  and  this  can  lead  to  some  wild 
oscillations  in  certain  circumstances  as  can  be  seen  in  the  example  in  Figure  20. 

The  last  technique  investigated  was  Akima  Cubic  Spline  Interpolation  (IMSL  routine  DCSAKM). 
This  routine  is  based  on  a  method  by  Akima  (1970)  to  combat  wiggles  in  the  interpolant  The 
result  is  that  the  shape  of  the  curve  produced  by  DCSAKM  matches  the  shape  of  the  data.  This, 
clearly,  is  what  is  required  and  produces  the  best  results  of  the  three  interpolation  techniques. 
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INTERPOLATION  OF  K(THETA) 


Linear  -  Natural  Cubic  Spline  -  Akima  Cubic  Spline 


Figure  20  :  Interpolations  of  K(0) 
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8.  Error  Messages 

The  following  is  a  list  of  possible  error  messages 

1100  Error  opening  parameter  file  "YVAKE.PAR" 

The  file  "WAKE.PAR"  must  be  in  the  same  directory  as  the  executable  file 

1101  Error  occurred  reading  one  of  the  parameters. 

This  means  that  an  error  occurred  reading  the  parameter  file,  possibly  one  of 
the  parameters  comment  lines  is  not  exactly  71  characters  or  the  comment 
line  contains  tabs 

1102  Invalid  NUMTHETA  in  Parameter  File  (0  <=  numtheta  <=  maxtheta) 
NUMTHETA  must  be  odd. 

The  number  of  theta  values  to  calculate  data  for  must  be  within  the  given 
range  and  must  be  odd. 

1103  Invalid  MODE  in  Parameter  File  (  0  <=  MODE  <=  maxmod). 

The  mode  to  calculate  data  for  must  be  within  the  given  range. 

1104  Invalid  THETAMAX  in  Parameter  File  (  THETAMAX  <  jt/2  ). 

The  maximum  9  value  allowed  must  be  less  than  7t/2.  A  value  of  1.56  is 
recommended. 

1105  Invalid  SUBDTH  in  Parameter  File  (0  <=  SUBDTH  <=  DEPTH  ). 

The  depth  of  the  submarine  must  be  within  the  given  range 

1106  Invalid  SUBSPD  in  Parameter  File  (0  <=  SUBSPD  <=  maxspd  ). 

The  speed  of  the  submarine  must  be  within  the  given  range 

1107  Invalid  SUBLEN  in  Parameter  File  (0  <=  SUBLEN  <=  msxien  ). 

The  length  of  the  submarine  must  be  within  the  given  range 

1108  Invalid  SUBRAD  in  Parameter  FUe  (0  <=  SUB  RAD  <=  maxrad  ). 

The  radius  of  the  submarine  must  be  within  the  given  range 
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1114 


1115 

1116 


1117 

1118 

1119 
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Invalid  XEND  in  Parameter  File  (  XEND  >=  XSTART  >=0  ) 

Invalid  XSTEP  in  Parameter  File  (  XSTEP  >  0  ). 

Invalid  YEND  in  Parameter  File  (  YEND  >=  YSTART  >=0  ) 

Invalid  YSTEP  in  Parameter  File  (  YSTEP  >  0  ). 

Invalid  ZSTART  in  Parameter  File  (  ZSTART  must  be  a  multiple  of 
"VALID  GRID  POINT" 

All  Z  values  must  lie  on  the  Z  Grid  as  defined  by  the  density  profile.  The  Z 
Grid  is  defined  by  "(((NUMDEN-l)/2)+l)"  Z  values  evenly  spaced  between 
0  and  DEPTH 

Invalid  ZEND  in  Parameter  File  (  ZEND  must  be  a  multiple  of 
"VALID  GRID  POINT" 

All  Z  values  must  lie  on  the  Z  Grid  as  defined  by  the  density  profile.  The  Z 
Grid  is  defined  by  "(((NUMDEN-l)/2)+l)"  Z  values  evenly  spaced  between 
0  and  DEPTH 

Invalid  ZEND  in  Parameter  File  (  ZEND  >=  ZSTART  >=0  ) 

Invalid  ZSTEP  in  Parameter  File  (  ZSTEP  must  be  a  multiple  of 
"VALID  GRID  POINT" 

All  Z  values  must  lie  on  the  Z  Grid  as  defined  by  the  density  profile.  The  Z 
Grid  is  defined  by  "(((NUMDEN-l)/2)+l)"  Z  values  evenly  spaced  between 
0  and  DEPTH 

Invalid  ZSTEP  in  Parameter  File  (  ZSTEP  >=0  ) 

Invalid  DEPTH  in  Parameter  File  (  DEPTH  >  0) 

The  depth  of  the  ocean  must  be  greater  than  zero 

Invalid  COMPU  in  Parameter  File  (COMPU  =  "Y"  or  "N"). 

Compute  U  Velocity  Component,  COMPU,  must  be  "Y"  for  Yes!,  or 
"N"  for  No! 

Invalid  COMPV  in  Parameter  File  (COMPV  =  "Y"  or  "N"). 


Compute  V  Velocity  Component,  COMPV,  must  be  "Y"  for  Yes!,  or 
"N"  for  No! 
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1121  Invalid  COMPW  in  Parameter  File  (COMPW  =  "Y"  or  ”N"). 

Compute  W  Velocity  Component,  COMPW,  must  be  "Y”  for  Yes!,  or 
"N"  for  No! 

1200  Error  occurred  opening  density  profile  file. 

Either  the  file  specified  by  the  parameter  DENFIL  is  missing  or  corrupt 

1201  Density  profile  file  has  a  non  odd  number  of  points. 

The  density  profile  must  have  an  odd  number  of  points 

1202  Density  Profile  contains  "X"  points,  not  "NUMDEN”  points  as  stated  in 
the  parameter  file. 

The  parameter  NUMDEN  in  the  Parameter  File  should  state  the  correct 
number  of  points  that  are  contained  in  the  density  profile.  DENFIL. 

1300  Error  opening  debugging  file. 

1400  Error  opening  velocity  component  file. 

1500  Error  opening  intermediate  integration  file 

1501  Error  occurred  writing  A(0)  data  to  the  integration  file 
Either  the  integration  file  is  missing  or  it  is  corrupt 

1700  Negative  Sigma  Value  passed  to  Function  LOOKUP. 

1701  Sigma  >  SigarrayCMAXMOD)  in  Function  LOOKUP. 

The  sigma  value  passed  to  the  function  lookup  is  greater  than  the  maximum 
sigma  value  allowed,  (i.e.  the  Sigma  value  for  K=0  for  the  maximum  mode 
calculated  (mode  30)). 
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